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Abstract 

The ability to discover physical laws and governing equations from data is one of hu¬ 
mankind's greatest intellectual achievements. A quantitative understanding of dynamic con¬ 
straints and balances in nature has facilitated rapid development of knowledge and enabled 
advanced technological achievements, including aircraft, combustion engines, satellites, and 
electrical power. In this work, we combine sparsity-promoting techniques and machine learn¬ 
ing with nonlinear dynamical systems to discover governing physical equations from measure¬ 
ment data. The only assumption about the structure of the model is that there are only a few 
important terms that govern the dynamics, so that the equations are sparse in the space of pos¬ 
sible functions; this assumption holds for many physical systems. In particular, we use sparse 
regression to determine the fewest terms in the dynamic governing equations required to accu¬ 
rately represent the data. The resulting models are parsimonious, balancing model complexity 
with descriptive ability while avoiding overfitting. We demonstrate the algorithm on a wide 
range of problems, from simple canonical systems, including linear and nonlinear oscillators 
and the chaotic Lorenz system, to the fluid vortex shedding behind an obstacle. The fluid ex¬ 
ample illustrates the ability of this method to discover the underlying dynamics of a system 
that took experts in the community nearly 30 years to resolve. We also show that this method 
generalizes to parameterized, time-varying, or externally forced systems. 

Keywords- Dynamical systems. Sparse regression. System identification. Compressed sensing. 


1 Introduction 

Extracting physical laws from data is a central challenge in many diverse areas of science and en¬ 
gineering. There are many critical data-driven problems, such as understanding cognition from 
neural recordings, inferring patterns in climate, determining stability of financial markets, predict¬ 
ing and suppressing the spread of disease, and controlling turbulence for greener transportation 
and energy. With abundant data and elusive laws, it is likely that data-driven discovery of dy¬ 
namics will continue to play an increasingly important role in these efforts. 

Advances in machine learning |[T9ll and data science 12^ |2T1 have promised a renaissance in 
the analysis and understanding of complex data, extracting patterns in vast multimodal data that 
is beyond the ability of humans to grasp. However, despite the rapid development of tools to 
understand static data based on statistical relationships, there has been slow progress in distill¬ 
ing physical models of dynamic processes from big data. This has limited the ability of data 
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science models to extrapolate the dynamics beyond the attractor where they were sampled and 
constructed. 

An analogy may be drawn with the discoveries of Kepler and Newton. Kepler, equipped with 
the most extensive and accurate planetary data of the era, developed a data-driven model for the 
motion of the planets, resulting in his famous elliptic orbits. However, this was an attractor based 
view of the world, and it did not explain the fundamental dynamic relationships that give rise 
to planetary orbits, or provide a model for how these bodies react when perturbed. Newton, in 
contrast, discovered a dynamic relationship between momentum and energy that described the 
underlying processes responsible for these elliptic orbits. This dynamic model may be generalized 
to predict behavior in regimes where no data was collected. Newton's model has proven remark¬ 
ably robust for engineering design, making it possible to land a spacecraft on the moon, which 
would not have been possible using Kepler's model alone. 

A seminal breakthrough by Schmidt and Lipson IHIlOl has resulted in a new approach to de¬ 
termine the underlying structure of a nonlinear dynamical system from data. This method uses 
symbolic regression (i.e., genetic programming |B2| ) to find nonlinear differential equations, and 
it balances complexity of the model, measured in the number of terms, with model accuracy. The 
resulting model identification realizes a long-sought goal of the physics and engineering commu¬ 
nities to discover dynamical systems from data. However, symbolic regression is expensive, does 
not scale well to large systems of interest, and may be prone to overfitting unless care is taken to 
explicitly balance model complexity with predictive power. In 1401 , the Pareto front is used to find 
parsimonious models in a large family of candidate models. 

In this work, we re-envision the dynamical system discovery problem from an entirely new 
perspective of sparse regression ll42l[T^[T8ll and compressed sensing IIT^l8l[9l[7ll2ll4^. In particular, 
we leverage the fact that most physical systems have only a few relevant terms that define the 
dynamics, making the governing equations sparse in a high-dimensional nonlinear function space. 
Before the advent of compressive sampling, and related sparsity-promoting methods, determining 
the few non-zero terms in a nonlinear dynamical system would have involved a combinatorial 
brute-force search, meaning that the methods would not scale to larger problems with Moore's 
law. However, powerful new theory guarantees that the sparse solution may be determined with 
high-probability using convex methods that do scale favorably with problem size. The resulting 
nonlinear model identification inherently balances model complexity (i.e., sparsity of right hand 
side dynamics) with accuracy, and the underlying convex optimization algorithms ensure that the 
method will be applicable to large-scale problems. 

The method described here shares some similarity to the recent dynamic mode decomposition 
(DMD), which is a linear dynamic regression 13611^ . DMD is an example of an equation-free 
method f 20ll , since it only relies on measurement data, but not on knowledge of the governing 
equations. Recent advances in the extended DMD have developed rigorous connections between 
DMD built on nonlinear observable functions and the Koopman operator theory for nonlinear 
dynamical systems lMil3Bll . However, there is currently no theory for which nonlinear observ¬ 
able functions to use, so that assumptions must be made on the form of the dynamical system. 
In contrast, the method developed here results in a sparse, nonlinear regression that automatically 
determines the relevant terms in the dynamical system. The trend to exploit sparsity in dynamical 
systems is recent but growing l38ll3^ 1251 1^ l35l[T]|. In this work, promoting sparsity in the dynam¬ 
ics results in parsimonious natural laws. 
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2 Background 

There is a long and fruitful history of modeling dynamics from data, resulting in powerful tech¬ 
niques for system identification 12^ . Many of these methods arose out of the need to understand 
complex flexible structures, such as the Hubble space telescope or the international space station. 
The resulting models have been widely applied in nearly every branch of engineering and applied 
mathematics, most notably for model-based feedback control. However, methods for system iden¬ 
tification typically require assumptions on the form of the model, and most often result in linear 
dynamics, limiting their effectiveness to small amplitude transient perturbations around a fixed 
point of the dynamics ITSl . 

This work diverges from the seminal work on system identification, and instead builds on 
symbolic regression and sparse representation. In particular, symbolic regression is used to find 
nonlinear functions describing the relationships between variables and measured dynamics (i.e., 
time derivatives). Traditionally, model complexity is balanced with describing capability using 
parsimony arguments such as the Pareto front. Here, we use sparse representation to determine 
the relevant model terms in an efficient and scalable framework. 

2.1 Symbolic regression and machine learning 

Symbolic regression involves the determination of a function that relates input-output data, and it 
may be viewed as a form of machine learning. Typically, the function is determined using genetic 
programming, which is an evolutionary algorithm that builds and tests candidate functions out of 
simple building blocks l22l . These functions are then modified according to a set of evolutionary 
rules and generations of functions are tested until a pre-determined accuracy is achieved. 

Recently, symbolic regression has been applied to data from dynamical systems, and ordinary 
differential equations were discovered from measurement data 1401 . Because it is possible to over¬ 
fit with symbolic regression and genetic programming, a parsimony constraint must be imposed, 
and in Il40l , they accept candidate equations that are at the Pareto front of complexity. 

2.2 Sparse representation and compressive sensing 

In many regression problems, only a few terms in the regression are important, and a sparse feature 
selection mechanism is required. For example, consider data measurements y G that may 
be a linear combination of columns from a feature library © G the linear combination of 

columns is given by entries of the vector ^ G so that: 

y = (1) 

Performing a standard regression to solve for ^ will result in a solution with nonzero contributions 
in each element. However, if sparsity of ^ is desired, so that most of the entries are zero, then it is 
possible to add an regularization term to the regression, resulting in the LASSO l[T4ilT^l42ll : 

I = argmin ||©|' - y ||2 + A|||'||i. (2) 

The parameter A weights the sparsity constraint. This formulation is closely related to the com¬ 
pressive sensing framework, which allows for the sparse vector ^ to be determined from rela¬ 
tively few incoherent random measurements IIT2l[8l[9l[7ll2ll4^. The sparse solution ^ to Eq.[^may 
also be used for sparse classification schemes, such as the sparse representation for classification 
(SRC) ||44|. Importantly, the compressive sensing and sparse representation architectures are con¬ 
vex and scale well to large problems, as opposed to brute-force combinatorial alternatives. 
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3 Sparse identification of nonlinear dynamics (SINDy) 


In this work, we are concerned with identifying the governing equations that underly a physical 
system based on data that may be realistically collected in simulations or experiments. Generically, 
we seek to represent the system as a nonlinear dynamical system 


±{t) = f(x(t)). 


(3) 


The vector x(t) = [a:i(t) X 2 {t) ■ ■ ■ Xn{t)]'^ E represents the state of the system at time t, 

and the nonlinear function f(x(t)) represents the dynamic constraints that define the equations of 
motion of the system. In the following sections, we will generalize Eq. ^ to allow the dynamics f 
to vary in time, and also with respect to a set of bifurcation parameters /x G M^. 

The key observation in this paper is that for many systems of interest, the function f often 
consists of only a few terms, making it sparse in the space of possible functions. For example, 
the Lorenz system in Eq. ( |22c| ) has very few terms in the space of polynomial functions. Recent 
advances in compressive sensing and sparse regression make this viewpoint of sparsity favorable, 
since it is now possible to determine which right hand side terms are non-zero without performing 
a computationally intractable brute-force search. 

To determine the form of the function f from data, we collect a time-history of the state x(t) 
and its derivative x(t) sampled at a number of instances in time ti, ^ 2 , • • • , tm- These data are then 
arranged into two large matrices: 


state 


rx^(ti)i 


Xi{ti) 

Xlit2) 

X2{tl) ■ 
X2{t2) ■ 

•• Xnih) 

■ ■ Xn{t2) 



_Xl {tm) 

X2itrn) 

Xnipm) _ 

ri^(ti)i 

i^(t2) 


Xl{tl) 

Xlit2) 

X2{tl) ■ ■ 

X2{t2) ■ • 

to 



Xl{tm) 

X2{tm) ■ • 



(4a) 


(4b) 


Next, we construct an augmented library ©(X) consisting of candidate nonlinear functions of the 
columns of X. For example, ©(X) may consist of constant, polynomial and trigonometric terms: 


0(X) 


1 X X^2 X^3 ••• sin(X) cos(X) sin(2X) cos(2X) 


(5) 


Here, higher polynomials are denoted as X^^ ^ ^Pz ^ For example, X^^ denotes the quadratic 

nonlinearities in the state variable x, given by: 


X^2 = 

xlih) 

xl{t2) 

Xl{tl)x2{tl) ■ 

Xl{t2)x2{t2) ■ 

•• xlih) 

•• xl(t2) 

X2{tl)xz{ti) ■ 

X2{t2)x?,{t2) ■ 

■■ xlih)' 
•• xlit2) 


_X\(im) 

X\(tm)X2{trn) 

X2(trn) 

X2itrn)Xz{tm) 
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Figure 1: Schematic of our algorithm for sparse identification of nonlinear dynamics, demonstrated on 
the Lorenz equations. Data is collected from measurements of the system, including a time history of the 
states X and derivatives X. Next, a library of nonlinear functions of the states, ©(X), is constructed. This 
nonlinear feature library is used to find the fewest terms needed to satisfy X = 0(X)S. The few entries 
in the vectors of S, solved for by sparse regression, denote the relevant terms in the right-hand side of the 
dynamics. Parameter values are a = 10,/d = 8/3, p = 28, (xq, Vo^zq)^ = (—8, 7,27)^. The trajectory on the 
Lorenz attractor is colored by the adaptive time-step required, with red requiring a smaller tilmestep. 


Each column of ©(X) represents a candidate function for the right hand side of Eq. There 
is tremendous freedom of choice in constructing the entries in this matrix of nonlinearities. Since 
we believe that only a few of these nonlinearities are active in each row of f, we may set up a 
sparse regression problem to determine the sparse vectors of coefficients H = ^2 ’ ’ ’ ^n] 

that determine which nonlinearities are active, as illustrated in Eig.[^ 

X = ©(X)H. (7) 

Each column of H represents a sparse vector of coefficients determining which terms are 
active in the right hand side for one of the row equations = 4(x) in Eq. Once H has been 
determined, a model of each row of the governing equations may be constructed as follows: 

Xfc = ffc(x) = 0(x'^)^fc. (8) 

Note that ©(x^) is a vector of symbolic functions of elements of x, as opposed to ©(X), which is 
a data matrix. This results in the overall model 

i = f(x) = S^(©(x^)f. (9) 

We may solve for H in Eq. Q using sparse regression. 
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3.1 Algorithm for sparse representation of dynamics with noise 

There are a number of algorithms to determine sparse solutions S to the regression problem in 
Eq. 0. Each column of Eq. 0 requires a distinct optimization problem to find the sparse vector 
of coefficients for the row equation. 

Eor the examples in this paper, the matrix ©(X) has dimensions m x p, where p is the number 
of candidate nonlinear functions, and where m ^ p since there are more time samples of data 
than there are candidate nonlinear functions. In most realistic cases, the data X and X will be 
contaminated with noise so that Eq. 0 does not hold exactly. In the case that X is relatively clean 
but the derivatives X are noisy, the equation becomes 

X = ©(X)H + 77Z, (10) 

where Z is a matrix of independent identically distributed Gaussian entries with zero mean, and 
p is the noise magnitude. Thus we seek a sparse solution to an overdetermined system with noise. 

The LASSO IHl [421 from statistics works well with this type of data, providing a sparse re¬ 
gression. However, it may be computationally expensive for very large data sets. 

An alternative is to implement the sequential thresholded least-squares algorithm in Code 0- 
In this algorithm, we start with a least-squares solution for H and then threshold all coefficients 
that are smaller than some cutoff value A. Once the indices of the remaining non-zero coefficients 
are identified, we obtain another least-squares solution for H onto the remaining indices. These 
new coefficients are again thresholded using A, and the procedure is continued until the non-zero 
coefficients converge. This algorithm is computationally efficient, and it rapidly converges to a 
sparse solution in a small number of iterations. The algorithm also benefits from simplicity, with 
a single parameter A required to determine the degree of sparsity in H. 

Depending on the noise, it may still be necessary to filter the data X and derivative X before 
solving for H. In particular, if only the data X is available, and X must be obtained by differentia¬ 
tion, then the resulting derivative matrix may have large noise magnitude. To counteract this, we 
use the total variation regularized derivative llQl to de-noise the derivative. An alternative would 
be to filter the data X and X, for example using the optimal hard threshold for singular values 
described in tT^ . 

It is important to note that previous algorithms to identify dynamics from data have been quite 
sensitive to noise HqI. The algorithm in Code 0 is remarkably robust to noise, and even when 
velocities must be approximated from noisy data, the algorithm works surprisingly well. 

Code 1: Sparse representation algorithm in Matlab. 

%% compute Sparse regression: sequential least squares 

Xi = Theta\dXdt ; % initial guess: Least-squares 

% lambda is our spars!float ion knob. 

for k=l:10 

smallinds = (abs (Xi ) <lambda) ; % find small coefficients 

Xi ( smallinds )=0; % and threshold 

for ind = l:n % n is state dimension 

biginds = ~ smallinds (:, ind) ; 

% Regress dynamics onto remaining terms to find sparse Xi 
Xi(biginds , ind) = Theta (:, biginds ) \dXdt (:, ind) ; 

end 

end 
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3.2 Cross-validation to determine parsimonious sparse solution on Pareto front 

To determine the sparsification parameter A in the algorithm in Code Q, it is helpful to use the 
concept of cross-validation from machine learning. It is always possible to hold back some test 
data apart from the training data to test the validity of models away from training values. In 
addition, it is important to consider the balance of model complexity (given by the number of 
nonzero coefficients in H) with the model accuracy. There is an "elbow" in the curve of accuracy 
vs. complexity parameterized by A, the so-called Pareto front. This value of A represents a good 
tradeoff between complexity and accuracy, and it is similar to the approach taken in 1401 . 

3.3 Extensions and Connections 

There are a number of extensions to the basic theory above that generalize this approach to a 
broader set of problems. First, the method is generalized to a discrete-time formulation, establish¬ 
ing a connection with the dynamic mode decomposition (DMD). Next, high-dimensional systems 
obtained from discretized partial differential equations are considered, extending the method to 
incorporate dimensionality reduction techniques to handle big data. Finally, the sparse regression 
framework is modified to include bifurcation parameters, time-dependence, and external forcing. 

3.3.1 Discrete-time representation 

The aforementioned strategy may also be implemented on discrete-time dynamical systems: 

XA; + 1 = f(XA;). (11) 

There are a number of reasons to implement Eq. ( pd] ). First, many systems, such as the logistic 
map in Eq. are inherently discrete-time systems. In addition, it may be possible to recover 
specific integration schemes used to advance Eq. The discrete-time formulation also foregoes 
the calculation of a derivative from noisy data. The data collection will now involve two matrices 
and X^: 


( 12 ) 


The continuous-time sparse regression problem in Eq. 0 now becomes: 

X^ = ©(X7^-^)H (13) 

and the function f is the same as in Eq. 0. 

In the discrete setting in Eq. ( [TT] >, and for linear dynamics, there is a striking resemblance to 
dynamic mode decomposition. In particular, if ©(x) = x, so that the dynamical system is linear, 
then Eq. ([0 becomes 

= ^ (X^)^ = 2^(X^-^)^. (14) 

This is equivalent to the DMD, which seeks a dynamic regression onto linear dynamics H^. In 
particular, is n x n dimensional, which may be prohibitively large for a high-dimensional state 
X. Thus, DMD identifies the dominant terms in the eigendecomposition of H^. 
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X2 

Xq 
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3.3.2 High-dimensional systems, partial differential equations, and dimensionality reduction 

Often, the physical system of interest may be naturally represented by a partial differential equa¬ 
tion (PDE) in a few spatial variables. If data is collected from a numerical discretization or from 
experimental measurements on a spatial grid, then the state dimension n may be prohibitively 
large. For example, in fluid dynamics, even simple two-dimensional and three-dimensional flows 
may require tens of thousands up to billions of variables to represent the discretized system. 

The method described above is prohibitive for a large state dimension n, both because of the 
factorial growth of © in n and because each of the n row equations in Eq. ^ requires a separate op¬ 
timization. Fortunately, many high-dimensional systems of interest evolve on a low-dimensional 
manifold or attractor that may be well-approximated using a dimensionally reduced low-rank ba¬ 
sis ^ csiiazi. For example, if data X is collected for a high-dimensional system as in Eq. ( [4a] ), it is 
possible to obtain a low-rank approximation using the singular value decomposition (SVD): 

= ^EV*. (15) 

In this case, the state x may be well approximated in a truncated modal basis given by the first 
r columns of ^ from the SVD: 


X ^ (16) 

where a is an r-dimensional vector of mode coefficients. We assume that this is a good approxi¬ 
mation for a relatively low rank r. Thus, instead of using the original high-dimensional state x, it 
is possible to obtain a sparse representation of the Galerkin projected dynamics fp in terms of the 
coefficients a: 

a = fp(a). (17) 

There are many choices for a low-rank basis, including proper orthogonal decomposition (POD) ||3l 
|l6l, based on the SVD. 

3.3.3 External forcing, bifurcation parameters, and normal forms 

In practice, many real-world systems depend on parameters, and dramatic changes, or bifurca¬ 
tions, may occur when the parameter is varied |[T5ll2^ . The algorithm above is readily extended 
to encompass these important parameterized systems, allowing for the discovery of normal forms 
associated with a bifurcation parameter /x. First, we append /x to the dynamics: 

X = f(x;/x) (18a) 

^ = 0. (18b) 

It is then possible to identify the right hand side f(x; /x) as a sparse combination of functions of 
components in x as well as the bifurcation parameter /x. This idea is illustrated on two examples, 
the one-dimensional logistic map and the two-dimensional Hopf normal form. 

Time-dependence, known external forcing or control u{t) may also be added to the dynamics: 

X = f (x, u{t ), t) (19a) 

i = 1. (19b) 

This generalization makes it possible to analyze systems that are externally forced or controlled. 
For example, the climate is both parameterized and has external forcing, including carbon 
dioxide and solar radiation. The financial market presents another important example with forc¬ 
ing and active feedback control, in the form of regulations, taxes, and interest rates. 
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4 Results 


We demonstrate the methods described in Sec. on a number of canonical systems, ranging from 
simple linear and nonlinear damped oscillators, to noisy measurements of the fully chaotic Lorenz 
system, and to measurements of the unsteady fluid wake behind a circular cylinder, extending this 
method to nonlinear partial differential equations (PDEs) and high-dimensional data. Finally, we 
show that bifurcation parameters may be included in the sparse models, recovering the correct 
normal forms from noisy measurements of the logistic map and the Hopf normal form. 

4.1 Example 1: Simple illustrative systems 

4.1.1 Example la: Two-dimensional damped oscillator (linear vs. nonlinear) 

In this example, we consider the two-dimensional damped harmonic oscillator with either linear 
or cubic dynamics, as in Eq. ( [^ . The dynamic data and the sparse identified model are shown in 
Fig. 1^ The correct form of the nonlinearity is obtained in each case; the augmented nonlinear li¬ 
brary ©(x) includes polynomials in x up to fifth order. The sparse identified model and algorithm 
parameters are shown in the Appendix in Tables and 



Figure 2: Comparison of linear system (left) and system with cubic nonlinearity (right). The sparse identi¬ 
fied system correctly identifies the form of the dynamics and accurately reproduces the phase portraits. 
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4.1.2 Example lb: Three-dimensional linear system 

A linear system with three variables and the sparse approximation are shown in Fig. In this 
case, the dynamics are given by 


X 


■- 0.1 -2 0 


X 

y 


2 - 0.1 0 


y 

z 


CO 

o 

1 

o 

o 

1 


z 


( 21 ) 


The sparse identification algorithm correctly identifies the system in the space of polynomials up 
to second or third order, and the sparse model is given in Table Interestingly, including poly¬ 
nomial terms of higher order (e.g. orders 4 or 5) introduces a degeneracy in the sparse identifica¬ 
tion algorithm, because linear combinations of powers of may approximate other exponential 
rates. This unexpected degeneracy motivates a hierarchical approach to identification, where sub¬ 
sequently higher order terms are included until the algorithm either converges or diverges. 




Figure 3: Three-dimensional linear system (solid colored lines) is well-captured by sparse identified system 
(dashed black line). 


4.2 Example 2: Lorenz system (Nonlinear ODE) 

Here, we consider the nonlinear Lorenz system fM l to explore the identification of chaotic dynam¬ 
ics evolving on an attractor, shown in Fig.[^ 


X 

= CF{y - x) 

(22a) 

y 

1 

TT 

1 

II 

(22b) 

z 

II 

1 

(22c) 


Although these equations give rise to rich and chaotic dynamics that evolve on an attractor, there 
are only a few terms in the right-hand side of the equations. Figure shows a schematic of how 
data is collected for this example, and how sparse dynamics are identified in a space of possible 
right-hand side functions using convex £i-minimzation. 

For this example, data is collected for the Lorenz system, and stacked into two large data 
matrices X and X, where each row of X is a snapshot of the state x in time, and each row of X 
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is a snapshot of the time derivative of the state x in time. Here, the right-hand side dynamics are 
identified in the space of polynomials 0(X) in (x, y, z) up to fifth order; 


©(X) 


x(t) y(t) z(t) x(t)2 x(t)y(t) x(t)z(t) y(t)^ 



(23) 


Each column of 0(X) represents a candidate function for the right hand side of Eq. and a 
sparse regression determines which terms are active in the dynamics, as in Fig.j^ and Eq. (0. 

Zero-mean Gaussian measurement noise with variance r] is added to the derivative calculation 
to investigate the effect of noisy derivatives. The short-time (t = 0 to t = 20) and long-time (t = 0 
to t = 250) system reconstruction is shown in Fig. for two different noise values, r] = 0.01 and 
7] = 10. The trajectories are also shown in dynamo view in Fig. and the £2 error vs. time for 
increasing noise r] is shown in Fig. Although the £2 error increases for large noise values 77 , 
the form of the equations, and hence the attractor dynamics, are accurately captured. Because 
the system has a positive Lyapunov exponent, small differences in model coefficients or initial 
conditions grow exponentially, until saturation, even though the attractor may remain unchanged. 

In the Lorenz example, the ability to capture dynamics on the attractor is more important than 
the ability to predict an individual trajectory, since chaos will quickly cause any small variations 
in initial conditions or model coefficients to diverge exponentially. As shown in Fig. our sparse 
model identification algorithm accurately reproduces the attractor dynamics from chaotic trajec¬ 
tory measurements. The algorithm not only identifies the correct linear and quadratic terms in the 
dynamics, but it accurately determines the coefficients to within .03% of the true values. When 
the derivative measurements are contaminated with noise, the correct dynamics are identified, 
and the attractor is well-preserved for surprisingly large noise values. When the noise is too large, 
the structure identification fails before the coefficients become too inaccurate. 


Full Simulation 





Figure 4: Trajectories of the Lorenz system for short-time integration from t = 0 to t = 20 (top) and 
long-time integration from t = Otot = 250 (bottom). The full dynamics (left) are compared with the sparse 
identified systems (middle, right) for various additive noise. The trajectories are colored by At, the adaptive 
Runge-Kutta time step. This color is a proxy for local sensitivity 
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For this example, we use the standard parameters a = 10,/3 = 8/3, p = 28, with an initial 
condition [x y = [—8 7 27]^. Data is collected from t = 0 to t = 100 with a time-step of 
At = 0.001. 

T ] — 0.01 77 = 10 






Figure 5: Dynamo view of trajectories of the Lorenz system. The exact system is shown in black (—) and 
the sparse identified system is shown in the dashed red arrow ( — ). 



Figure 6: Error vs. time for sparse identified systems generated from data with increasing sensor noise y. 
This error corresponds to the difference between solid black and dashed red curves in Fig. Sensor noise 
values are y G {0.0001,0.001,0.01,0.1,1.0,10.0}. 
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4.3 Example 3: Fluid wake behind a cylinder (Nonlinear PDF) 

The Lorenz system is a low-dimensional model of more realistic high-dimensional partial differ¬ 
ential equation (PDF) models for fluid convection in the atmosphere. Many systems of interest are 
governed by PDFs IIMI . such as weather and climate, epidemiology, and the power grid, to name 
a few. Fach of these examples are characterized by big data, consisting of large spatially resolved 
measurements consisting of millions or billions of states and spanning orders of magnitude of 
scale in both space and time. However, many high-dimensional, real-world systems evolve on a 
low-dimensional attractor, making the effective dimension much smaller fT^ . 

Here we generalize the sparse identification of nonlinear dynamics method to an example 
in fluid dynamics that typifies many of the challenges outlined above. Data is collected for the 
fluid flow past a cylinder at Reynolds number 100 using direct numerical simulations of the two- 
dimensional Navier-Stokes equations ilHUTTII . Then, the nonlinear dynamic relationship between 
the dominant coherent structures is identified from these flow field measurements with no knowl¬ 
edge of the governing equations. 

The low-Reynolds number flow past a cylinder is a particularly interesting example because 
of its rich history in fluid mechanics and dynamical systems. It has long been theorized that 
turbulence may be the result of a sequence of Hopf bifurcations that occur as the Reynolds number 
of the flow increases 1^ . The Reynolds number is a rough measure of the ratio of inertial and 
viscous forces, and an increasing Reynolds number may correspond, for example, to increasing 
flow velocity, giving rise to more rich and intricate structures in the fluid. 

After 15 years, the first Hopf bifurcation was discovered in a fluid system, in the transition 
from a steady laminar wake to laminar periodic vortex shedding at Reynolds number 47 IITtI 
Il5l|32l. This discovery led to a long-standing debate about how a Hopf bifurcation, with cubic 
nonlinearity, can be exhibited in a Navier-Stokes fluid with quadratic nonlinearities. After 15 
more years, this was finally resolved using a separation of time-scales argument and a mean-field 
model Il3n , shown in Fq. ( [24| ). It was shown that coupling between oscillatory modes and the 
base flow gives rise to a slow manifold (see Fig. left), which results in algebraic terms that 
approximate cubic nonlinearities on slow timescales. 



A - vortex shedding 


- POD mode 1 





B - mean flow 


Uy - POD mode 2 




C - unstable fixed point 

Uz - shift mode 




’few 


(feler - t; 




1 0 1 2 3 4 5 6 7! 
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Figure 7: Illustration of the low-rank dynamics underlying the periodic vortex shedding behind a circular 
cylinder at low Reynolds number. Re = 100. 
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This example provides a compelling test-case for the proposed algorithm, since the under¬ 
lying form of the dynamics took nearly three decades to uncover. Indeed, the sparse dynamics 
algorithm correctly identifies the on-attractor and off-attractor dynamics using quadratic nonlin¬ 
earities and preserves the correct slow-manifold dynamics. It is interesting to note that when the 
off-attractor trajectories are not included in the system identification, the algorithm incorrectly 
identifies the dynamics using cubic nonlinearities, and fails to correctly identify the dynamics 
associated with the shift mode, which connects the mean flow to the unstable steady state. 

4.3.1 Direct numerical simulation 

The direct numerical simulation involves a fast multi-domain immersed boundary projection 
method lUTl [TTII . Four grids are used, each with a resolution of 450 x 200, with the finest grid 
having dimensions of 9 x 4 cylinder diameters and the largest grid having dimensions of 72 x 32 
diameters. The finest grid has 90,000 points, and each subsequent coarser grid has 67,500 distinct 
points. Thus, if the state includes the vorticity at each grid point, then the state dimension is 
292,500. The vorticity field on the finest grid is shown in Fig.[^ The code is non-dimensionalized 
so that the cylinder diameter and free-stream velocity are both equal to one: D = 1 and Uoo = h 
respectively. The simulation time-step is At = 0.02 non dimensional time units. 


4.3.2 Mean field model 

To develop a mean-field model for the cylinder wake, first we must reduce the dimension of 
the system. The proper orthogonal decomposition (POD) ||T6|, provides a low-rank basis that is 
optimal in the sense, resulting in a hierarchy of orthonormal modes that are ordered by mode 
energy. The first two most energetic POD modes capture a significant portion of the energy; the 
steady-state vortex shedding is a limit cycle in these coordinates. An additional mode, called the 
shift mode, is included to capture the transient dynamics connecting the unstable steady state 
with the mean of the limit cycle ll^ (i.e., the direction connecting point 'C' to point 'B' in Fig.[^. 

In the three-dimensional coordinate system described above, the mean-field model for the 
cylinder dynamics are given by: 


X = fix — ujy + Axz 
y = fiy Ayz 

z = —X{z — x‘^ — y‘^) 


(24a) 

(24b) 

(24c) 


If A is large, so that the z-dynamics are fast, then the mean flow rapidly corrects to be on the (slow) 
manifold z = x^^ + y‘^ given by the amplitude of vortex shedding. When substituting this algebraic 
relationship into Eqs. |24a and |24b[ we recover the Hopf normal form on the slow manifold. 

Remarkably, similar dynamics are discovered by the sparse dynamics algorithm, purely from 
data collected from simulations. The identified model coefficients, shown in Table only include 
quadratic nonlinearities, consistent with the Navier-Stokes equations. Moreover, the transient 
behavior, shown in Figs. and [TO} is captured qualitatively for solutions that do not start on the 
slow manifold. When the off-attractor dynamics in Fig. are not included in the training data, the 
model incorrectly identifies a simple Hopf normal form in x and y with cubic nonlinearities. 

The data from Fig. was not included in the training data, and although qualitatively similar, 
the identified model does not exactly reproduce the transients. Since this initial condition had 
twice the fluctuation energy in the x and y directions, the slow manifold approximation may not 
be valid here. Relaxing the sparsity condition, it is possible to obtain models that agree almost 
perfectly with the data in Figs.[8p0} although the model includes higher order nonlinearities. 
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Full Simulation 


Identified System 



200 



Figure 8: Evolution of the cylinder wake trajectory in reduced coordinates. The full simulation (left) comes 
from direct numerical simulation of the Navier-Stokes equations, and the identified system (right) captures 
the dynamics on the slow manifold. Color indicates simulation time. 



Figure 9: Evolution of the cylinder wake trajectory starting from a flow state initialized at the mean of 
the steady-state limit cycle. Both the full simulation and sparse model capture the off-attractor dynamics, 
characterized by rapid attraction of the trajectory onto the slow manifold. 


Full Simulation 



Identified System 



Figure 10: This simulation corresponds to an initial condition obtained by doubling the magnitude of the 
limit cycle behavior. This data was not included in the training of the sparse model. 
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4.4 Example 4: Bifurcations and Normal Forms 

It is also possible to identify normal forms associated with a bifurcation parameter fi by suspend¬ 
ing it in the dynamics as a variable; 

X = f(x;/i) (25a) 

A = 0. (25b) 

It is then possible to identify the right hand side f (x; p) as a sparse combination of functions of 

components in x as well as the bifurcation parameter p. This idea is illustrated on two examples, 
the one-dimensional logistic map and the two-dimensional Hopf normal form. 

4.4.1 Logistic map 

The logistic map is a classical model that exhibits a cascade of bifurcations, leading to chaotic 
trajectories. The dynamics with stochastic forcing tjk and parameter n are given by 

Xk+i = fiXkil - Xk) + Vk- (26) 

Sampling the stochastic system at ten parameter values of /i, the algorithm correctly identifies the 
underlying parameterized dynamics, shown in Fig.pTjand Table 





Figure 11: Attracting sets of the logistic map vs. the parameter ja. (left) Data from stochastically 
forced system and (right) the sparse identified system. Data is sampled at rows indicated in red for 
/i G {2.5,2.75,3,3.25,3.5,3.75,3.8,3.85,3.9,3.95}. The forcing r]k is Gaussian with magnitude 0.025. 
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4.4.2 Hopf normal form 


The final example illustrating the ability of the sparse dynamics method to identify parameterized 
normal forms is the Hopf normal form 1281 . Noisy data is collected from the Hopf system 

X = fix + ujy — Axix^^ + y^) (27a) 

y = -ujx + yy- Ay{x^ + y^) (27b) 


for various values of the parameter y. Data is collected on the blue and red trajectories in Fig. 12 


and noise is added to simulate sensor noise. The total variation derivative ITOll is used to de-noise 
the derivative for use in the algorithm. 

The sparse model identification algorithm correctly identifies the Hopf normal form, with 
model parameters given in Table The noise-free model reconstruction is shown in Fig.[^ Note 
that with noise in the training data, although the model terms are correctly identified, the actual 
values of the cubic terms are off by almost 8%. Collecting more training data or reducing the noise 
magnitude both improve the model agreement. 



Figure 12: Training data to identify Hopf normal form. Blue trajectories denote solutions that start outside 
of the fixed point for /i < 0 or the limit cycle for /i > 0, and red trajectories denote solutions that start inside 
of the limit cycle. 
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5 Discussion 


In summary, we have demonstrated a powerful new technique to identify nonlinear dynamical 
systems from data without assumptions on the form of the governing equations. This builds on 
prior work in symbolic regression but with innovations related to sparse regression, which al¬ 
low our algorithms to scale to high-dimensional systems. We demonstrate this new method on 
a number of example systems exhibiting chaos, high-dimensional data with low-rank coherence, 
and parameterized dynamics. As shown in the Lorenz example, the ability to predict a specific 
trajectory may be less important than the ability to capture the attractor dynamics. The example 
from fluid dynamics highlights the remarkable ability of this method to extract dynamics in a 
fluid system that took three decades for experts in the community to explain. There are numerous 
fields where this method may be applied, where there is ample data and the absence of governing 
equations, including neuroscience, climate science, epidemiology, and financial markets. Fields 
that already use genetic programming, such as machine learning control for turbulent fluid sys¬ 
tems 01311/ rnay also benefit. Finally, normal forms may be discovered by including parameters 
in the optimization, as shown on two examples. The identification of sparse governing equations 
and parameterizations marks a significant step toward the long-held goal of intelligent, unassisted 
identification of dynamical systems. 

A number of open problems remain surrounding the dynamical systems aspects of this pro¬ 
cedure. For example, many systems possess dynamical symmetries and conserved quantities that 
may alter the form of the identified dynamics. For example, the degenerate identification of a 
linear system in a space of high-order polynomial nonlinearities suggest a connection with near¬ 
identity transformations and dynamic similarity. We believe that this may be a fruitful line of 
research. Finally, it will be important to identify which approximating function space to use based 
on the data available. For example, it may be possible to improve the function space to make the 
dynamics more sparse through subsequent coordinate transformations IT^ . 

Data science is not a panacea for all problems in science and engineering, but used in the 
right way, it provides a principled approach to maximally leverage the data that we have and 
inform what new data to collect. Big data is happening all across the sciences, where the data 
is inherently dynamic, and where traditional approaches are prone to overfitting. Data discovery 
algorithms that produce parsimonious models are both rare and desirable. Data-science will only 
become more critical to efforts in science in engineering, where data is abundant, but physical 
laws remain elusive. These efforts include understanding the neural basis of cognition, extracting 
and predicting coherent changes in the climate, stabilizing financial markets, managing the spread 
of disease, and controlling turbulence. 
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Appendix 


Table 1; Damped harmonic oscillator with linear terms. 


r 

'xdot' 

'ydot' 

1 ' 

[ 0] 

[ 0] 

x' 

[-0.1015] 

[-1.9990] 

Y' 

[ 2.0027] 

[-0.0994] 

XX ' 

[ 0] 

[ 0] 

xy' 

[ 0] 

[ 0] 

yy' 

[ 0] 

[ 0] 

xxx' 

[ 0] 

[ 0] 

xxy' 

[ 0] 

[ 0] 

xyy' 

[ 0] 

[ 0] 

yyy' 

[ 0] 

[ 0] 

xxxx' 

[ 0] 

[ 0] 

xxxy' 

[ 0] 

[ 0] 

xxyy' 

[ 0] 

[ 0] 

xyyy' 

[ 0] 

[ 0] 

yyyy' 

[ 0] 

[ 0] 

xxxxx' 

[ 0] 

[ 0] 

xxxxy' 

[ 0] 

[ 0] 

xxxyy' 

[ 0] 

[ 0] 

xxyyy' 

[ 0] 

[ 0] 

xyyyy' 

[ 0] 

[ 0] 

yyyyy' 

[ 0] 

[ 0] 
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Table 2: Damped harmonic oscillator with cubic nonlinearity. 


r 

'xdot' 

'ydot' 

1 ' 

[ 0] 

[ 0] 

x' 

[ 0] 

[ 0] 

Y' 

[ 0] 

[ 0] 

XX ' 

[ 0] 

[ 0] 

xy' 

[ 0] 

[ 0] 

yy' 

[ 0] 

[ 0] 

xxx' 

[-0.0996] 

[-1.9994] 

xxy' 

[ 0] 

[ 0] 

xyy' 

[ 0] 

[ 0] 

yyy' 

[ 1.9970] 

[-0.0979] 

xxxx' 

[ 0] 

[ 0] 

xxxy' 

[ 0] 

[ 0] 

xxyy' 

[ 0] 

[ 0] 

xyyy' 

[ 0] 

[ 0] 

yyyy' 

[ 0] 

[ 0] 

xxxxx' 

[ 0] 

[ 0] 

xxxxy' 

[ 0] 

[ 0] 

xxxyy' 

[ 0] 

[ 0] 

xxyyy' 

[ 0] 

[ 0] 

xyyyy' 

[ 0] 

[ 0] 

yyyyy' 

[ 0] 

[ 0] 


Table 3: Three-dimensional linear system. 


r 

' xdot' 

'ydot' 

'zdot' 

1 ' 

[ 0] 

[ 0] 

[ 0] 

x' 

[-0.0996] 

[-1.9997] 

[ 0] 

Y' 

[ 2.0005] 

[-0.0994] 

[ 0] 

z' 

[ 0] 

[ 0] 

[-0.3003] 

XX ' 

[ 0] 

[ 0] 

[ 0] 

xy' 

[ 0] 

[ 0] 

[ 0] 

xz' 

[ 0] 

[ 0] 

[ 0] 

yy' 

[ 0] 

[ 0] 

[ 0] 

yz' 

[ 0] 

[ 0] 

[ 0] 

zz' 

[ 0] 

[ 0] 

[ 0] 
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Table 4; Lorenz system identified using sparse representation with rj = 1.0. 


r 

'xdot' 

' ydot' 

'zdot' 

1' 

[ 0] 

[ 0] 

[ 0] 

x' 

[-9.9996] 

[27.9980] 

[ 0] 

Y' 

[ 9.9998] 

[-0.9997] 

[ 0] 

z' 

[ 0] 

[ 0] 

[-2.6665] 

XX ' 

[ 0] 

[ 0] 

[ 0] 

xy' 

[ 0] 

[ 0] 

[ 1.0000] 

xz' 

[ 0] 

[-0.9999] 

[ 0] 

yy' 

[ 0] 

[ 0] 

[ 0] 

yz' 

[ 0] 

[ 0] 

[ 0] 

zz' 

[ 0] 

[ 0] 

[ 0] 

xxx' 

[ 0] 

[ 0] 

[ 0] 

xxy' 

[ 0] 

[ 0] 

[ 0] 

xxz' 

[ 0] 

[ 0] 

[ 0] 

xyy' 

[ 0] 

[ 0] 

[ 0] 

xyz' 

[ 0] 

[ 0] 

[ 0] 

xzz' 

[ 0] 

[ 0] 

[ 0] 

yyy' 

[ 0] 

[ 0] 

[ 0] 

yyz' 

[ 0] 

[ 0] 

[ 0] 

yzz' 

[ 0] 

[ 0] 

[ 0] 

zzz' 

[ 0] 

[ 0] 

[ 0] 

xxxx' 

[ 0] 

[ 0] 

[ 0] 

xxxy' 

[ 0] 

[ 0] 

[ 0] 

xxxz' 

[ 0] 

[ 0] 

[ 0] 

xxyy' 

[ 0] 

[ 0] 

[ 0] 

xxyz' 

[ 0] 

[ 0] 

[ 0] 

xxzz' 

[ 0] 

[ 0] 

[ 0] 

xyyy' 

[ 0] 

[ 0] 

[ 0] 

xyyz' 

[ 0] 

[ 0] 

[ 0] 

xyzz' 

[ 0] 

[ 0] 

[ 0] 

xzzz' 

[ 0] 

[ 0] 

[ 0] 

yyyy' 

[ 0] 

[ 0] 

[ 0] 

yyyz' 

[ 0] 

[ 0] 

[ 0] 

yyzz' 

[ 0] 

[ 0] 

[ 0] 

yzzz' 

[ 0] 

[ 0] 

[ 0] 

zzzz' 

[ 0] 

[ 0] 

[ 0] 

xxxxx' 

[ 0] 

[ 0] 

[ 0] 

xxxxy' 

[ 0] 

[ 0] 

[ 0] 

xxxxz' 

[ 0] 

[ 0] 

[ 0] 

xxxyy' 

[ 0] 

[ 0] 

[ 0] 

xxxyz' 

[ 0] 

[ 0] 

[ 0] 

xxxzz' 

[ 0] 

[ 0] 

[ 0] 

xxyyy' 

[ 0] 

[ 0] 

[ 0] 

xxyyz' 

[ 0] 

[ 0] 

[ 0] 

xxyzz' 

[ 0] 

[ 0] 

[ 0] 

xxzzz' 

[ 0] 

[ 0] 

[ 0] 

xyyyy' 

[ 0] 

[ 0] 

[ 0] 

xyyyz' 

[ 0] 

[ 0] 

[ 0] 

xyyzz' 

[ 0] 

[ 0] 

[ 0] 

xyzzz' 

[ 0] 

[ 0] 

[ 0] 

xzzzz' 

[ 0] 

[ 0] 

[ 0] 

yyyyy' 

[ 0] 

[ 0] 

[ 0] 

YYYYz' 

[ 0] 

[ 0] 

[ 0] 

yyyzz' 

[ 0] 

[ 0] 

[ 0] 

yyzzz' 

[ 0] 

[ 0] 

[ 0] 

yzzzz' 

[ 0] 

[ 0] 

[ 0] 

zzzzz' 

[ 0] 

[ 0] 

[ 0] 


23 



Table 5: Dynamics of cylinder wake modes using sparse representation. Notice that quadratic terms are 
identified. 


r 1 

' xdot' 

'ydot' 


' zdot' 


1' 

[ -0.1225] 

[ -0. 

. 0569] 

[ -20.8461] 

x' 

[ -0.0092] 

[ 1. 

. 0347] 

[-4.6476e- 

-04] 

Y' 

[ -1.0224] 

[ 0. 

.0047] 

[ 2.4057e- 

-04] 

z' 

[-9.2203e-04] 

[-4.4932e-04] 

[ -0.2968] 

XX ' 

[ 0] 

[ 

0] 

[ 0.0011] 

xy' 

[ 0] 

[ 

0] 

[ 

0] 

xz' 

[ 2.1261e-04] 

[ 0. 

.0022] 

[ 

0] 

yy' 

[ 0] 

[ 

0] 

[ 8.6432e- 

-04] 

yz' 

[ -0.0019] 

[ -0. 

.0018] 

[ 

0] 

zz' 

[ 0] 

[ 

0] 

[ -0.0010] 

xxx' 

[ 0] 

[ 

0] 

[ 

0] 

xxy' 

[ 0] 

[ 

0] 

[ 

0] 

xxz' 

[ 0] 

[ 

0] 

[ 

0] 

xyy' 

[ 0] 

[ 

0] 

[ 

0] 

xyz' 

[ 0] 

[ 

0] 

[ 

0] 

xzz' 

[ 0] 

[ 

0] 

[ 

0] 

yyy' 

[ 0] 

[ 

0] 

[ 

0] 

yyz' 

[ 0] 

[ 

0] 

[ 

0] 

yzz' 

[ 0] 

[ 

0] 

[ 

0] 

zzz' 

[ 0] 

[ 

0] 

[ 

0] 

xxxx' 

[ 0] 

[ 

0] 

[ 

0] 

xxxy' 

[ 0] 

[ 

0] 

[ 

0] 

xxxz' 

[ 0] 

[ 

0] 

[ 

0] 

xxyy' 

[ 0] 

[ 

0] 

[ 

0] 

xxyz' 

[ 0] 

[ 

0] 

[ 

0] 

xxzz' 

[ 0] 

[ 

0] 

[ 

0] 

xyyy' 

[ 0] 

[ 

0] 

[ 

0] 

xyyz' 

[ 0] 

[ 

0] 

[ 

0] 

xyzz' 

[ 0] 

[ 

0] 

[ 

0] 

xzzz' 

[ 0] 

[ 

0] 

[ 

0] 

yyyy' 

[ 0] 

[ 

0] 

[ 

0] 

yyyz' 

[ 0] 

[ 

0] 

[ 

0] 

yyzz' 

[ 0] 

[ 

0] 

[ 

0] 

yzzz' 

[ 0] 

[ 

0] 

[ 

0] 

zzzz' 

[ 0] 

[ 

0] 

[ 

0] 

xxxxx' 

[ 0] 

[ 

0] 

[ 

0] 

xxxxy' 

[ 0] 

[ 

0] 

[ 

0] 

xxxxz' 

[ 0] 

[ 

0] 

[ 

0] 

xxxyy' 

[ 0] 

[ 

0] 

[ 

0] 

xxxyz' 

[ 0] 

[ 

0] 

[ 

0] 

xxxzz' 

[ 0] 

[ 

0] 

[ 

0] 

xxyyy' 

[ 0] 

[ 

0] 

[ 

0] 

xxyyz' 

[ 0] 

[ 

0] 

[ 

0] 

xxyzz' 

[ 0] 

[ 

0] 

[ 

0] 

xxzzz' 

[ 0] 

[ 

0] 

[ 

0] 

xyyyy' 

[ 0] 

[ 

0] 

[ 

0] 

xyyyz' 

[ 0] 

[ 

0] 

[ 

0] 

xyyzz' 

[ 0] 

[ 

0] 

[ 

0] 

xyzzz' 

[ 0] 

[ 

0] 

[ 

0] 

xzzzz' 

[ 0] 

[ 

0] 

[ 

0] 

yyyyy' 

[ 0] 

[ 

0] 

[ 

0] 

YYYYz' 

[ 0] 

[ 

0] 

[ 

0] 

yyyzz' 

[ 0] 

[ 

0] 

[ 

0] 

yyzzz' 

[ 0] 

[ 

0] 

[ 

0] 

yzzzz' 

[ 0] 

[ 

0] 

[ 

0] 

zzzzz' 

[ 0] 

[ 

0] 

[ 

0] 
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Table 6: Logistic map identified using sparse representation. 


r 

'x_{k+l}' 

'r_{k+l} 

1' 

[ 0] 

[ 0] 

x' 

[ 0] 

[ 0] 

r' 

[ 0] 

[1.0000] 

XX ' 

[ 0] 

[ 0] 

xr' 

[ 0.9993] 

[ 0] 

rr' 

[ 0] 

[ 0] 

xxx' 

[ 0] 

[ 0] 

xxr' 

[-0.9989] 

[ 0] 

xrr' 

[ 0] 

[ 0] 

rrr' 

[ 0] 

[ 0] 

xxxx' 

[ 0] 

[ 0] 

xxxr' 

[ 0] 

[ 0] 

xxrr' 

[ 0] 

[ 0] 

xrrr' 

[ 0] 

[ 0] 

rrrr' 

[ 0] 

[ 0] 

xxxxx' 

[ 0] 

[ 0] 

xxxxr' 

[ 0] 

[ 0] 

xxxrr' 

[ 0] 

[ 0] 

xxrrr' 

[ 0] 

[ 0] 

xrrrr' 

[ 0] 

[ 0] 

rrrrr' 

[ 0] 

[ 0] 
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Table 7: Hopf normal form identified using sparse representation. Here u represents the bifurcation pa¬ 
rameter /i. 


r 

'xdot' 

' ydot' 

'udot' 

1 ' 

[ 0] 

[ 0] 

[ 0] 

x' 

[ 0] 

[ 0.9914] 

[ 0] 

Y' 

[-0.9920] 

[ 0] 

[ 0] 

u' 

[ 0] 

[ 0] 

[ 0] 

XX ' 

[ 0] 

[ 0] 

[ 0] 

xy' 

[ 0] 

[ 0] 

[ 0] 

xu' 

[ 0.9269] 

[ 0] 

[ 0] 

yy' 

[ 0] 

[ 0] 

[ 0] 

yu' 

[ 0] 

[ 0.9294] 

[ 0] 

uu' 

[ 0] 

[ 0] 

[ 0] 

xxx' 

[-0.9208] 

[ 0] 

[ 0] 

xxy' 

[ 0] 

[-0.9244] 

[ 0] 

xxu' 

[ 0] 

[ 0] 

[ 0] 

xyy' 

[-0.9211] 

[ 0] 

[ 0] 

xyu' 

[ 0] 

[ 0] 

[ 0] 

xuu' 

[ 0] 

[ 0] 

[ 0] 

yyy' 

[ 0] 

[-0.9252] 

[ 0] 

yyu' 

[ 0] 

[ 0] 

[ 0] 

yuu' 

[ 0] 

[ 0] 

[ 0] 

uuu' 

[ 0] 

[ 0] 

[ 0] 

xxxx' 

[ 0] 

[ 0] 

[ 0] 

xxxy' 

[ 0] 

[ 0] 

[ 0] 

xxxu' 

[ 0] 

[ 0] 

[ 0] 

xxyy' 

[ 0] 

[ 0] 

[ 0] 

xxyu' 

[ 0] 

[ 0] 

[ 0] 

xxuu' 

[ 0] 

[ 0] 

[ 0] 

xyyy' 

[ 0] 

[ 0] 

[ 0] 

xyyu' 

[ 0] 

[ 0] 

[ 0] 

xyuu' 

[ 0] 

[ 0] 

[ 0] 

xuuu' 

[ 0] 

[ 0] 

[ 0] 

yyyy' 

[ 0] 

[ 0] 

[ 0] 

yyyu' 

[ 0] 

[ 0] 

[ 0] 

yyuu' 

[ 0] 

[ 0] 

[ 0] 

yuuu' 

[ 0] 

[ 0] 

[ 0] 

uuuu' 

[ 0] 

[ 0] 

[ 0] 

xxxxx' 

[ 0] 

[ 0] 

[ 0] 

xxxxy ' 

[ 0] 

[ 0] 

[ 0] 

xxxxu' 

[ 0] 

[ 0] 

[ 0] 

xxxyy ' 

[ 0] 

[ 0] 

[ 0] 

xxxyu' 

[ 0] 

[ 0] 

[ 0] 

xxxuu' 

[ 0] 

[ 0] 

[ 0] 

xxyyy' 

[ 0] 

[ 0] 

[ 0] 

xxyyu' 

[ 0] 

[ 0] 

[ 0] 

xxyuu' 

[ 0] 

[ 0] 

[ 0] 

xxuuu' 

[ 0] 

[ 0] 

[ 0] 

xyyyy' 

[ 0] 

[ 0] 

[ 0] 

xyyyu' 

[ 0] 

[ 0] 

[ 0] 

xyyuu' 

[ 0] 

[ 0] 

[ 0] 

xyuuu' 

[ 0] 

[ 0] 

[ 0] 

xuuuu' 

[ 0] 

[ 0] 

[ 0] 

yyyyy' 

[ 0] 

[ 0] 

[ 0] 

YYYYU' 

[ 0] 

[ 0] 

[ 0] 

yyyuu' 

[ 0] 

[ 0] 

[ 0] 

yyuuu' 

[ 0] 

[ 0] 

[ 0] 

yuuuu' 

[ 0] 

[ 0] 

[ 0] 

uuuuu' 

[ 0] 

[ 0] 

[ 0] 
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